Emergence of phenotypic and genotypic antimicrobial resistance in Mycobacterium tuberculosis

Concentration dependency of phenotypic and genotypic isoniazid-rifampicin resistance emergence was investigated to obtain a mechanistic understanding on how anti-mycobacterial drugs facilitate the emergence of bacterial populations that survive throughout treatment. Using static kill curve experiments, observing two evolution cycles, it was demonstrated that rifampicin resistance was the result of non-specific mechanisms and not associated with accumulation of drug resistance encoding SNPs. Whereas, part of isoniazid resistance could be accounted for by accumulation of specific SNPs, which was concentration dependent. Using a Hollow Fibre Infection Model it was demonstrated that emergence of resistance did not occur at concentration–time profiles mimicking the granuloma. This study showed that disentangling and quantifying concentration dependent emergence of resistance provides an improved rational for drug and dose selection although further work to understand the underlying mechanisms is needed to improve the drug development pipeline.

www.nature.com/scientificreports/ wall permeability. Antibiotic tolerance refers to an increased ability of the bacterial population to survive antibiotic exposure without increase in MIC. These manifests itself in greater minimum duration to kill 99% of the population (MDK 99% ). Antibiotic tolerance is caused by a combination of general and drug specific mechanisms that are independent of antibiotic class and include reduced growth rates, metabolic shifts and increased activity of efflux pumps 6,7 . Antibiotic persistence refers to a predestined sub-population of bacteria, often 0.01-1% of an inherently heterogenous population of bacteria, that is able to survive throughout antibiotic treatment. The persistent sub-population is more difficult to kill, i.e. has a longer MDK 99% , while the majority of the bacterial population is fully susceptible, i.e. has a shorter MDK 99% , and this results in a characteristic bi-phasic killing profile over the course of antibiotic exposure without increase in MIC. Antibiotic persistence is not inheritable, meaning a heterogenous population will regrow when antibiotic pressure is taken off [6][7][8] . Bacteria within the second phase of a bi-phasic killing profile may consist of different sub-populations. Little research has been done in the domain of disentangling and quantifying genotypic from phenotypic resistance in the tail of bi-phasic killing profiles for antimicrobial combination therapy, even though it forms the cornerstone of anti-TB therapy. Isoniazid and rifampicin form the backbone of standard drug-sensitive anti-tuberculosis therapy and isoniazid enhances bi-phasic bacillary killing from sputum, but concomitant therapy with rifampicin does not prevent this from happening 8,9 . Given the mutation rate, in H37Rv, of isoniazid (2.56 × 10 -8 -3.2 × 10 -710, 11 ) and rifampicin to a lesser extent (6 × 10 -10 -2.4 × 10 -711-15 ), this could be emergence of heterogenous genotypic resistance, phenotypic resistance or a combination of both. Disentangling and quantifying emergence of genotypic and phenotypic resistance over the course of isoniazid and rifampicin combination therapy can therefore inform ongoing in-vivo investigations evaluating the increase in rifampicin and isoniazid dose in order to achieve improved and sustained overall anti-mycobacterial activity 4 .
To that end the concentration dependency of emergence of phenotypic and genotypic isoniazid-rifampicin resistance was investigated and the impact of in-vivo mimicking pharmacokinetic profiles at a range of dose scenarios were explored.

Results
Pharmacodynamic interactions and emergence of resistance. Static kill curve experiments, to elucidate emergence of resistance defined by phenotype and genotype, comprised three stages; during the first evolution cycle, cultures were exposed to various isoniazid and/or rifampicin concentrations (0.25, 1 and 8 MIC equivalents) for one week, subsequently cultures were regrown in drug free media for three weeks and then entered to a second evolution cycle in which cultures were exposed to the same drug conditions as during the first evolution cycle for another week. This allowed for selection of bacterial sub-populations that were able to survive isoniazid and/or rifampicin exposure. Isoniazid and rifampicin drug effects disappear over the course of the experiment, whether given alone or in combination, both in the first and in the second evolution cycle.
However, when isoniazid was used alone, potency decreased (associated with an increase in IC 50 , i.e. the concentration at which half-maximum antimycobacterial effect is achieved, of 402%) after repeated exposure but a resistant population, i.e. described by the tail in bi-phasic mycobacterial killing model, emerged later in time compared to the first exposure cycle. This was associated with a decrease in model parameter λ, (λ is defined as onset of susceptibility loss as a function of time) of 16.9% (Table 1, Fig. 1A,C and Supplementary Fig. S1). Likewise, when studied alone rifampicin potency decreased (associated with an increase in IC 50 of 182%) after www.nature.com/scientificreports/ repeated exposure but unlike the isoniazid treatment resistance emerged around the same time in the first and second evolution cycle (Table 1 and Fig. 1B,D). After repeated concomitant exposure to rifampicin, susceptibility to isoniazid dropped through lower efficacy (associated with a decrease in E MAX , the maximum antimycobacterial effect, of 56.8%) ( Table 1 and Fig. 1A). Isoniazid resistance emerged later (associated with a decrease in model parameter λ of 32.0%) in the presence of rifampicin (Table 1 and Fig. 1C), i.e. the tail in bi-phasic mycobacterial killing curve emerged later. Presence of isoniazid improved rifampicin susceptibility (associated with an increase in E MAX of 42.1%) but rifampicin resistance emerged earlier (associated with an increase in model parameter λ of 58.2%) in the presence of isoniazid (Table 1 and Fig. 1B,D).
A bi-phasic killing profile was observed, that was more pronounced during the second evolution cycle, and this could consist of different genotypic resistant and phenotypic persistent sub-populations.
Emergence of genotypic resistance. Whole Genome Sequencing (WGS) of samples from the static kill curve experiments at baseline and at the end of the first and second evolution cycle was used to disentangle emergence of genotypic and phenotypic resistance. Emergence of phenotypic isoniazid resistance (Table 1 and Fig. 1C) coincides with appearance of SNPs in genes associated with isoniazid resistance, katG, inhA, fabG1, and kasA and these accumulate after repeated exposure ( Fig. 2A). However, katG isoniazid resistance-associated nucleotide changes occurred stochastically in isoniazid containing experiments while accumulation after repeated exposure was displayed (Fig. 2C). The frequency of isoniazid resistance-associated nucleotide changes was higher in experiments with lower isoniazid exposure and presence of rifampicin did not result in suppression of isoniazid resistance-associated nucleotide changes from emerging (Fig. 2C). This indicates that isoniazid induced bi-phasic mycobacterial killing can be partly explained by concentration dependent emergence of isoniazid resistant SNPs while the remainder can be attributed to emergence of phenotypic resistance.
No rifampicin resistance-associated nucleotide changes were observed and the mutation rate in genes associated with rifampicin resistance, i.e. rpoB and rpoC, was lower when compared to mutation rates in genes associated with isoniazid resistance ( Fig. 2A,B). Consequently, rifampicin induced bi-phasic mycobacterial killing could only be attributed to emergence of phenotypic resistance.
Impact of dosing regimen. The impact of various dosing strategies was evaluated using the Hollow-Fibre Infection Model (HFIM), an experimental setup that allows in-vivo like concentration-time profiles at the site of infection to be mimicked in-vitro 16   www.nature.com/scientificreports/ MGIT-TTP/day. However, there was no trend of bi-phasic bacillary clearance for any of the four tested dosing regimens using a HFIM (Fig. 3A) and drug-resistant associated SNPs only occurred at the start in one out of four experiments (TDS at 1/3 standard dose) and disappeared over the course of the experiment (Fig. 3B). This indicates that emergence of genotypic resistance did not occur at the mimicked dosing regimens.

Discussion
While investigating pharmacodynamic drug-drug interactions and emergence of resistance, with static kill curve experiments describing two evolution cycles, it was shown that rifampicin resistance was solely phenotypic (i.e. non-specific resistance mechanisms) but isoniazid could be attributed to phenotypic (non-specific) as well as specific genotypic mechanisms (Figs. 1, 2). Genotypic isoniazid resistance tends to occur at lower exposures (Fig. 2). However, for each of the four tested dosing regimens mimicked in the HFIM, that ranged from standard regimens 2 to increased dosing regimens similar to those currently being investigated in clinical trials 4,5 , there were no sign of emergence of resistance (Fig. 3). Since our findings also displayed antimicrobial exposure dependent bacillary clearance, as observed in patients 4,17 , we confirm that pharmacokinetic endpoints such as AUC/MIC and Cmax/MIC, i.e. the main pharmacokinetic endpoints in anti-tuberculosis therapy 18,19 , have been attained in all four dosing regimens and that exposures are well above target levels or emergence of resistance. However, previous studies with isoniazid and rifampicin in the HFIM did report the emergence of at least phenotypic resistance in the first 7 days [20][21][22] . A possible explanation for discrepancies in emergence of phenotypic resistance might be quantification methods for bacillary load given that inoculums were similar around 10 6 . Previous studies used Colony Forming Units (CFU's) as the quantification method, while here we used MGIT-TTP, which is known to have increased sensitivity associated with optimized liquid growth media and detection of oxygen consumption 23 . Nevertheless, our results show that emergence of genotypic resistance at clinical and investigational dose levels is unlikely provided patients adhere to their treatment regimen. However, more detailed investigations could provide an improved mechanistic understanding of the benefits to be expected from optimization of dosing regimens 4,5 , and accounting for genetic polymorphism in genes coding for metabolizing enzymes 24 . This research would require focus on gene expression, metabolomics, proteomics and lipidomics in order to dissect and quantify sub-populations by shifts in lipid metabolism, cell wall thickening as well as drug specific responses such as mycolic acid pathway caused by isoniazid and upregulation of drug targets for rifampicin 7 .
Adherence scenarios were not studied using the HFIM but on the basis of results from the static kill curve experiments it could be concluded that continued exposure to sub-or MIC equivalent drug levels increase the chance of genotypic isoniazid resistance emergence. This emphasises the importance of adherence to the treatment 25 in order to prevent low isoniazid exposure and consequent emergence of resistance. In this study only two evolution cycles were performed; a larger number of evolution cycles might have also resulted in emergence of rifampicin resistance as the mutation rate to rifampicin resistance is lower [11][12][13][14][15] compared to isoniazid 10,11 .
In this study we adopted the HFIM to model the dynamic exposure of bacteria to antimicrobials. There are discrepancies between static kill curve and HFIM models and these are likely a result of the experimental design, for example the persistence of the antibiotic in the media. It was confirmed that isoniazid and rifampicin crossed the fibers and not stick to plastics in the HFIM (Supplementary Fig. S2), differences are therefore probably a result of static experiments being susceptible to chemical degradation of antibiotics and depletion of nutrients. Isoniazid is chemically unstable at 37 °C and concentrations decrease by over 50% over 7 days 26 . While in the www.nature.com/scientificreports/ HFIM a new bolus was supplied every 24 or 8 h, for O.D. and T.D.S dosing respectively, isoniazid was not topped up over the course of a 7-day static kill curve evolution cycle. Emergence of phenotypic and genotypic isoniazid resistance in the static kill curve experiments might therefore also have been exacerbated by low isoniazid levels towards the end of an evolution cycle. Rifampicin has been shown to decrease by over 90% over 7 days at 37 °C and likewise phenotypic resistance emerged in static kill curve experiments might also have been facilitated in part by declining rifampicin levels over the course of the experiment. Furthermore, the volumes of the experimental models might have been a confounding factor, the volume of a tube in static kill curve experiments was only 8 mL and not refreshed over the course of a 7-day evolution cycle. While the system volume of the HFIM was 108 mL and refreshed on a continuous basis. The phenotypic resistance that emerged in rifampicin static kill curve experiments might have been stress and not drug associated which could explain the absence of rifampicin genotypic resistance.
Here, we demonstrate that emergence of antimicrobial specific (genotypic) resistance only occurs when antibiotic levels fall below MIC levels, however, MICs are maintained following clinical dosing provided that adherence to the regimen is good. Further work to understand the mechanisms underlying concentration dependent emergence of phenotypic resistance is desirable in order to improve the drug development pipeline and design of novel regimens. For example, dose reduction of linezolid when given with bedaquiline and pretomanid to treat drug resistant TB as considered recently by the WHO 27 .

Materials and methods
All static kill curve and hollow-fibre infection model experiments were conducted in biosafety level 3 laboratories.

Static kill-curve experiments.
All drug experiments and growth controls were performed in duplicate at 37 °C and ambient air. Isoniazid (40 mg/mL) and rifampicin (83 mg/L) stock solutions were prepared with sterile distilled water directly from the drug vials (Becton Dickinson-MGIT 960 SIRE KIT). Stock solutions were then diluted to drug experiment conditions containing 250, 1000 or 8000 ng/mL of isoniazid and 12.5, 50, and 400 ng/mL rifampicin alone or in combination with 10% OADC-supplemented Middlebrook broth 7H9 (Becton Dickinson-BBL MGIT 7ML).
M. tuberculosis H37Rv (NCTC 7416, ATCC 9360, obtained from Public Health England culture collections), was incubated in a MGIT Tubes containing 7 mL Middlebrook broth 7H9 (Becton Dickinson-BBL MGIT 7ML), supplemented with 10% OADC (Becton Dickinson-MGIT OADC ENRICHMENT 6 VIALS) prior to the first evolution cycle. There was a master tube for each experiment and evolution cycle one was started once the bacterial load of the master culture reached 10 5.5 CFU/mL. Each drug experiment or growth control experiment started with inoculating MGIT tubes containing 7 mL Middlebrook broth 7H9 with pellets from the master tube resuspended in 0.5 mL 7H9, 0.8 mL OADC and 0.2 mL drug or blank solution, rendering a total volume of 8.5 mL with a bacterial load of ~ 10 5.5 CFU/mL.
The bacteria were separated from media by centrifugation at 2683 Relative centrifugal force (RCF) for 10 min at the end of the first evolution cycle and the regrowth phase of the experiment. Pellets were washed in Sterile phosphate-buffered saline (PBS) and pelleted at 2683 RCF for 10 min twice. Thereafter the pellets were resuspended in the new condition at the start of the second or third stage in volumes as described for initial inocula.
Daily samples of 50 μL were taken for bacterial load quantification and three samples (baseline, after the first and after the second evolution cycle) were taken for Whole Genome Sequencing. www.nature.com/scientificreports/

Hollow-fibre infection model experiments.
A cellulosic (C3008, FibreCellSystems Ltd) cartridge was inoculated with 20 mL 10 5.5 CFU/mL M. tuberculosis, strain H37Rv. Incubation followed the identical procedure as per kill-curve experiments. A drug-free incubation phase of 2-3 days preceded the start of the drug experiments. Isoniazid and rifampicin concentration-time profiles mimicked unbound concentrations in the granuloma under the assumption that there is equilibrium between unbound drug concentrations in plasma and the granuloma 28 . Model predictions were obtained with a previously published model, adjusted for 42% and 83% plasma protein binding for isoniazid and rifampicin, respectively 29 . The following pharmacokinetic properties were mimicked for the OD and TDS standard dosage experiments: isoniazid C max = 3.32 mg/l and rifampicin C max = 0.788 mg/l with t 1/2 = 4.7 h). For the TDS experiment with 1/3 of standard dosing isoniazid C max was 1.06 mg/l and rifampicin C max = 0.263 mg/l with a t 1/2 at 4.7 h. The OD experiment at 3 times standard dosing simulated isoniazid C max at 9.50 mg/l and rifampicin C max at 2.36 mg/l with a t 1/2 at 4.7 h. A web application (https:// pkpdia. shiny apps. io/ hfs_ app/) 28 was subsequently used to convert secondary pharmacokinetic parameter estimates, CMAX/C0 and t1/2, into pump settings at a system volume of 108 mL (central reservoir, 50 mL; intracapillary space and tubing, 44 mL; and extracapillary space, 14 mL). First-order absorption of the drugs, to mimic oral absorption characteristics, was omitted and replaced by bolus administration in the central reservoir in these experiments to avoid complex experimental setups.
Growth control and drug experiments lasted for 14 and 7 days, respectively. Drug concentrations in the hollow-fibre medium were not measured during the experiments but pump-settings were validated using bacteria free experiments to ensure in-vivo mimicking pharmacokinetic profiles were simulated in-vitro ( Supplementary  Fig. S2).
Daily samples of 50 μL were taken for bacterial load quantification and two samples, at baseline at the end of the drug experiments, were taken for Whole Genome Sequencing.
Bacterial load quantification. Fifty μL samples from kill-curve and hollow-fibre infection model experiments were inoculated in fresh MGIT tubes containing 7 mL Middlebrook broth 7H9 with 0.8 mL OADC. MGIT-TTP was used as measure of bacterial load 23 .
Whole genome sequencing. The CTAB method was used to extract genomic DNA and Qubit dsDNA kits (Life Technologies) were used to quantify DNA prior to sequencing 30 . WGS for start culture and evolution cycle two samples were performed using an Illumina HiSeq platform and samples from the first evolution cycle and all Hollow-fibre infection model experiments were performed using the Illumina NextSeq platform. WGS was performed following the manufacturers' instructions and local validated protocols and results have been deposited (Supplementary Table S1).  Fig. S3 and Supplementary Data file S1). Model estimates were computed using the SAEM estimation method in nlmixr 2.0.4 on a Windows 10 operating system. Minus twice the log likelihood of the data was used as objective function value (OFV) and a drop in OFV of at least 3.84 (P = 0.05) was considered to improve the model's ability to fit the data with statistical significance after inclusion of one degree of freedom to a nested hierarchical model. Assessment of model performance was further supported by goodness-of-fit diagnostics including observation-population predictions, observations-individual predictions, Normalised Prediction Distribution Errors (NPDE)-population predictions and NPDE-time.
Baseline MGIT-TTP at experiment level ( P i ) was estimated using a typical baseline MGIT-TTP ( θ TV ) and a deviation from the typical baseline MGIT-TTP described by between experiment variability (η) (Eq. 1).
Data from the growth control experiments were used to develop a log-growth model (Eq. 2) describing changes in MGIT TTP bacterial load (B) over time using the parameters net growth ( θ k net ) and maximum carrying capacity ( θ B MAX ). The parameters θ k net and θ B MAX followed log-normal distribution for between experiment variability as described in Eq. (1). Data from re-growth experiments was subsequently included and predictive performance of the models was evaluated.
Isoniazid and rifampicin anti-mycobacterial effects were parameterised (Eq. 3) using maximum drug efficacy ( θ E MAX ) and potency ( θ IC 50 ), i.e. the concentration at which half-maximum inhibition is achieved, with C being the isoniazid or rifampicin concentration and ( θ γ ) being the shape factor.
(2) www.nature.com/scientificreports/ Loss of susceptibility to isoniazid and or rifampicin over the course of the experiment was parameterised using θ BETA , i.e. magnitude of susceptibility loss, and θ , i.e. onset of susceptibility loss as a function of time (Eq. 4). D r u g c o m b i n at i o n s w e r e e v a l u at e d a s a d d i t i v e d r u g e f f e c t s at f i r s t , i . e . Effect = (Effect isoniazid × Effect time isoniazid ) + (Effect rifampicin × Effect time rifampicin )) , and subsequently the Bliss independence model was tested and retained 31 , i.e. Effect = ((Effect isoniazid × Effect time isoniazid ) + (Effect rifampicin ×Effect time rifampicin )) − ((Effect isoniazid × Effect time isoniazid ) × (Effect rifampicin × Effect time rifampicin )) . Synergism and antagonism were evaluated as categorical effects for presence of drugs in combination and a drop in OFV of at least 3.84 (P = 0.05) was considered to improve the model's ability to fit the data with statistical significance in stepwise covariate model building 32 . Whole genome sequence analysis. Analysis of WGS data was divided in two steps. First, all mutations relative to the H37Rv reference genome were determined for 33 genes of interest associated with resistance to rifampicin, isoniazid, ethambutol, pyrazinamide, streptomycin, ethionamide, amikacin, capreomycin, kanamycin, fluoroquinolones, para-aminosalicylic acid, cycloserine, linezolid, bedaquiline, clofazimine and delamanid (rpoB, rpoC, fabG1, inhA, katG, kasA, ahpC, embR, embC, embA, embB, rpsA, pncA, panD, rpsL, gid, rrs, ethA, ethR, tlyA, eis, gyrB, gyrA, folC, ribD, thyX, thyA, ald, alr, rplC, rrl, Rv0678 and fbiA). Bcftools (v. 1.11), bwa (v. 0.7.12) and samtools (v. 1.9) were used for variant calling using the number of forward reference, reverse reference, forward non-reference, reverse non-reference alleles. The percentage of base pairs changed within genes, were stratified by genes associated with isoniazid and/or rifampicin resistance (katG, inhA, fabG1, kasA, rpoB and rpoC). Only base pairs with ≥ 15% changed to the H37Rv reference genome were used to calculate the percentage of base pairs changed within genes.
Tb-profiler (v. 3.0.1) was used to determine the SNPs by nucleotide changes 33 . Results were stratified by drug resistant nucleotide changes, i.e. isoniazid and rifampicin associated resistance, and other resistance that associate with any other of the aforementioned drugs.

Data availability
Accession numbers to WGS data are on an ENA project (PRJEB53040 & Supplementary Table S1) and the pharmacodynamic drug-drug interaction model is attached as Supplementary Data file S1. (3) Effect = θ E MAX × C θ γ θ IC 50 θ γ + C θ γ .